NB: this Rmd can be used for all types of ratings: emotion, arousal, valence
Abbreviations: - sub : subjects(s) - rat : rating(s)
library(tidyverse)
library(future)
library(furrr)
library(tictoc)
library(RNifti)
library(proxy) # distances
library(profvis)
library(DT)
library(formattable)
library(janitor)
library(ComplexHeatmap)
library(circlize) # For color mapping functions
library(grid) # For gpar()
# library(pheatmap)
# library(heatmaply)
library(papayaWidget)
library(ggstatsplot)
source("funs_V12_ROIs.R")
metadata_copes <- tribble(
~cope, ~emotion,
1, "anger",
2, "disgust",
3, "fear",
4, "happy",
5, "neutral",
6, "pain",
7, "sad"
)
# Do you want to remove neutral cope(s) from the RSA analysis?
# If "YES", the results filename will have option _RM (remove neutral)
remove_neutrals = "NO"
# choose the type of betas to be used following this correspondence table:
# rsa_ROI : one_ev_per_movie[_minus_neutral]
# rsa_emotion_predictors_ROI : emotion_predictors[_minus_neutral]
# rsa_emotion_high_low_predictors_ROI : emotion_high_low_predictors[_minus_neutral]
copes_type <- "emotion_predictors_minus_neutral"
# Choose one set of ROIs - e.g. an atlas - from ${bd_atlases} (ROIS_REPO)
atlas_filename <- "Yeo17.nii.gz"
# Schaefer.nii.gz
# Schaefer.nii.gz
# Yeo17.nii.gz
# Yeo17.nii.gz
# Atlases repos:
# https://www.lead-dbs.org/helpsupport/knowledge-base/atlasesresources/cortical-atlas-parcellations-mni-space/
# https://github.com/neurodata/neuroparc
# https://www.fmrib.ox.ac.uk/datasets/brainmap+rsns/
bd="/data00/leonardo/RSA/analyses"
bd_ratings = paste0(bd,"/RATINGS")
bd_atlases = paste0(bd,"/ROIS_REPO")
bd_results = paste0(bd,"/RSA_ROI_APP/results_RSA_ROI")
ratings_type <- c("emotion","arousal","valence","aroval")
rsa_flavour="rsa_emotion_predictors_ROI"
# Choose the distance metric to use for fmri RDM
# The metric for ratings RDM *should be* euclidean, since arousal and valence
# ratings have only one value
# Supported methods: 'pearson', 'spearman', 'euclidean', 'cosine' or 'mahalanobis'
dist_method_rating <- "euclidean"
dist_method_fmri = "euclidean"
dist_method_rsa = "pearson"
# Vector of zeropadded sub_ids
subs_file <- "/data00/leonardo/RSA/sub_list.txt"
subs <- sprintf("%02d",readLines(subs_file) %>% as.numeric)
# ---------- ONLY TOP RATERS BELOW ------------
subs <- c("02","03","12","11","22","09","29","28","26","32","23","15","20","19")
Extract the pathname of all copes using the list.files()
function. Also define a copes_numba vector with all the copes
numbers.
NB: The cope numbers in the cope column are NOT
zeropadded since this is how they come out from FSL Feat
# load df_path_copes and add metadata (label and emotion/neutral)
df_path_copes <- import_df_path_copes(bd, copes_type, rsa_flavour) %>%
inner_join(metadata_copes, by="cope") %>%
relocate(sub, cope, emotion, path)
# If remove_neutrals is "YES", remove neutral copes and arrange by sub/cope (just to be sure)
if (remove_neutrals == "YES") {
df_path_copes <- df_path_copes %>%
filter(emotion != "neutral") %>%
arrange(sub,cope)
}
# ---------- SPECIFIC TO EMOTION PREDICTION -------------
# add a label colum to the df_path_copes which represents the emotion
# (for the original 56 movies it was high_low_code, e.g. JK_A_high)
df_path_copes$label <- df_path_copes$emotion
# ---------- SPECIFIC TO EMOTION PREDICTION -------------
copes_numba <- df_path_copes$cope %>% unique
cat(
"there are",length(copes_numba)," copes including ",
nrow(df_path_copes %>% filter(sub==subs[1], emotion=="neutral"))," neutral copes \n"
)
## there are 7 copes including 1 neutral copes
# df_path_copes
# length(copes_numba)
df_path_copes
## # A tibble: 182 × 5
## sub cope emotion path label
## <chr> <dbl> <chr> <chr> <chr>
## 1 02 1 anger /data00/leonardo/RSA/analyses/emotion_predictors_m… anger
## 2 02 2 disgust /data00/leonardo/RSA/analyses/emotion_predictors_m… disg…
## 3 02 3 fear /data00/leonardo/RSA/analyses/emotion_predictors_m… fear
## 4 02 4 happy /data00/leonardo/RSA/analyses/emotion_predictors_m… happy
## 5 02 5 neutral /data00/leonardo/RSA/analyses/emotion_predictors_m… neut…
## 6 02 6 pain /data00/leonardo/RSA/analyses/emotion_predictors_m… pain
## 7 02 7 sad /data00/leonardo/RSA/analyses/emotion_predictors_m… sad
## 8 03 1 anger /data00/leonardo/RSA/analyses/emotion_predictors_m… anger
## 9 03 2 disgust /data00/leonardo/RSA/analyses/emotion_predictors_m… disg…
## 10 03 3 fear /data00/leonardo/RSA/analyses/emotion_predictors_m… fear
## # ℹ 172 more rows
atlas_path <- paste0(bd_atlases,"/",atlas_filename)
atlas_nii <- readNifti(paste0(bd_atlases,"/",atlas_filename))
region_labels <- atlas_nii[atlas_nii > 0] %>% unique %>% sort
# NB: ratings_type is now a *vector*, e.g.:
# ratings_type <- c("emotion","arousal","valence","aroval")
# to test the fn below for one rating type : do_RDMs_rats("emotion")
do_RDMs_rats <- function(ratings_type, remove_neutrals) {
ratings_path <- paste0(bd_ratings,"/",ratings_type,"_ratings.csv")
rats <- read_csv(ratings_path)
# ------- IMPORTANT : AVG ACROSS MOVIES FOR EMOTION PREDICTORS MODEL -----
rats <- rats %>%
select(sub, emotion, starts_with("r_")) %>%
group_by(sub, emotion) %>%
reframe(across(starts_with("r_"), median, na.rm = TRUE)) %>%
ungroup %>%
arrange(sub, emotion)
# ------- IMPORTANT : AVG ACROSS MOVIES FOR EMOTION PREDICTORS MODEL -----
if (remove_neutrals == "YES") {
rats <- rats %>% filter(emotion != "neutral")
}
rdm <- rats %>%
filter(sub %in% subs) %>%
select(sub, starts_with("r_")) %>%
group_by(sub) %>%
nest %>%
mutate(!!paste0("RDM_",ratings_type) := data %>% map(~ DDOS_vec(.x, dist_method_rating))) %>%
ungroup %>%
# remove the sub with the prospect of merging rdms of different ratings/models
select(!data)
return(rdm)
}
# calculate rdms for all ratings_type and return one column for each ratings_type
RDMs_rats <- map_dfc(ratings_type, ~ do_RDMs_rats(.x, remove_neutrals)) %>%
# remove the duplicated sub_ columns (generated by do_RDMs_rats)
janitor::clean_names() %>%
mutate(sub = sub_1) %>%
select(sub, starts_with("rdm"))
## Warning: There was 1 warning in `reframe()`.
## ℹ In argument: `across(starts_with("r_"), median, na.rm = TRUE)`.
## ℹ In group 1: `sub = "02"` and `emotion = "anger"`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))
# # example heatmap
# plot_tril(
# RDMs_rats[1,]$rdm_emotion, model = "emotion",
# fontsize=10, side_mm=80, reord = FALSE
# )
IMPORTANT : Since we already removed neutral from df_path_copes we don’t need to filter anything anymore here.
# Fn to calculate the rdm_fmri for one sub
# Test for one sub: do_RDM_fmri(subs[1], copes_numba, df_path_copes)
do_RDM_fmri <- function(sub_id, copes_numba, df_path_copes) {
# load all the copes nii.gz into a df
df_copes <- load_sub_copes(sub_id, copes_numba, df_path_copes)
# the code below does the calculation of the tril of RDM_fmri for all rois for one sub
one_sub_RDM_fmri <- tibble(
sub = sub_id,
roi = region_labels
) %>%
# extract idx_roi for atlas voxels in that roi
mutate(idx_roi = roi %>% map(~ which(atlas_nii == .x)) ) %>%
# extract the df_copes for the voxels in that region (idx_roi)
mutate(df_copes_region = idx_roi %>% map(~ df_copes[.x,]) ) %>%
# calculate tril of RDM_roi (output as numeric vector)
mutate(rdm_fmri = df_copes_region %>% map(~ DDOS_vec(t(.x), dist_method_fmri)) ) %>%
# select only sub, roi (numba) and RDM_roi
select(sub, roi, rdm_fmri)
return(one_sub_RDM_fmri)
}
# Calculate RDMs_fmri for all subs
plan(multisession, workers = 5)
RDMs_fmri <- subs %>% future_map_dfr(
~ {
paste0("Calculating RDMs for sub ",.x,"\n") %>% cat
do_RDM_fmri(.x, copes_numba, df_path_copes)
}
)
## Calculating RDMs for sub 02
## Calculating RDMs for sub 03
## Calculating RDMs for sub 12
## Calculating RDMs for sub 11
## Calculating RDMs for sub 22
## Calculating RDMs for sub 09
## Calculating RDMs for sub 29
## Calculating RDMs for sub 28
## Calculating RDMs for sub 26
## Calculating RDMs for sub 32
## Calculating RDMs for sub 23
## Calculating RDMs for sub 15
## Calculating RDMs for sub 20
## Calculating RDMs for sub 19
plan(sequential)
# RDMs_fmri
# The dist_method_rsa is set above.
# can be pearson, spearman
# dist_method_rsa <- "spearman"
# Join fmri and rats RDMs and create a copy of the rat RDM for each sub/roi
RSA <- right_join(RDMs_fmri, RDMs_rats, by = "sub") %>%
group_by(sub) %>%
# Dynamically create rsa_ columns for every ratings type
mutate(across(
all_of(paste0("rdm_", ratings_type)),
~ map2_dbl(.x, rdm_fmri, ~ cor(.x, .y, method = dist_method_rsa)),
.names = "rsa_{.col}"
)) %>%
mutate(across(starts_with("rsa_"), round, 3)) %>%
ungroup()
## Warning: There were 17 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `across(...)`.
## ℹ In group 6: `sub = "15"`.
## Caused by warning in `cor()`:
## ! the standard deviation is zero
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 16 remaining warnings.
# Get the mean RSA for each ratings_type
RSA_mean <- RSA %>%
select(roi, starts_with("rsa_")) %>%
group_by(roi) %>%
reframe(
across(
starts_with("rsa_"),
list(mean = ~ round(mean(.x, na.rm = TRUE), 2))
)
) %>%
rename_with(
~ str_replace(.x, "rsa_rdm_", ""),
starts_with("rsa_")
)
RSA_mean %>% datatable()
# Takes a df and a "column" containing a nested vector in each row
# and returns the element-wise mean of that vector across rows
get_elementwise_mean <- function(mydf, nested_col, fun=mean) {
m <- matrix(
mydf[[nested_col]] %>% unlist,
nrow = nrow(mydf),
ncol = length(mydf[[nested_col]][[1]]),
byrow = TRUE
)
m_summary <- apply(m, 2, fun)
return(m_summary)
}
# mean RDMs_rats
mean_RDMs_rats <- RDMs_rats %>% select(starts_with("rdm_")) %>% colnames %>%
map_dfc(~{
mean_tril <- get_elementwise_mean(RDMs_rats, .x, fun = mean)
tibble(!!.x := mean_tril)
})
# mean RDMs_fmri
plan(multisession, workers = 5)
mean_RDMs_fmri <- RDMs_fmri$roi %>%
future_map_dfr(~{
mean_tril <- get_elementwise_mean(RDMs_fmri %>% filter(roi == .x), "rdm_fmri")
tibble(roi = .x, mean_tril = mean_tril)
}) %>%
group_by(roi) %>%
nest
plan(sequential)
# Plot mean RDMs_rats
ratings_type %>%
map(
~ plot_tril(
mean_RDMs_rats[[paste0("rdm_",.x)]], model = .x,
fontsize=7, side_mm=80, reord = FALSE
)
)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
mean_RDMs_fmri
## # A tibble: 17 × 2
## # Groups: roi [17]
## roi data
## <dbl> <list>
## 1 1 <tibble [294 × 1]>
## 2 2 <tibble [294 × 1]>
## 3 3 <tibble [294 × 1]>
## 4 4 <tibble [294 × 1]>
## 5 5 <tibble [294 × 1]>
## 6 6 <tibble [294 × 1]>
## 7 7 <tibble [294 × 1]>
## 8 8 <tibble [294 × 1]>
## 9 9 <tibble [294 × 1]>
## 10 10 <tibble [294 × 1]>
## 11 11 <tibble [294 × 1]>
## 12 12 <tibble [294 × 1]>
## 13 13 <tibble [294 × 1]>
## 14 14 <tibble [294 × 1]>
## 15 15 <tibble [294 × 1]>
## 16 16 <tibble [294 × 1]>
## 17 17 <tibble [294 × 1]>
seq(nrow(mean_RDMs_fmri)) %>%
map(
~ plot_tril(
mean_RDMs_fmri[.x,]$data, model = paste0("fmri_roi_",.x),
fontsize=7, side_mm=80, reord = FALSE
)
)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
##
## [[15]]
##
## [[16]]
##
## [[17]]
RSA_mean %>% datatable()
roi_numba = 1
RSA %>%
select(sub, roi, starts_with("rsa_")) %>%
filter(roi == roi_numba) %>%
select(starts_with("rsa_")) %>%
pivot_longer(cols = starts_with("rsa_"), names_to = "model") %>%
mutate(model = str_replace(model, "rsa_rdm_","")) %>%
ggwithinstats(
x = model, y = value,
type = "p", # p, np, r, b
results.subtitle = TRUE
)
prepare_papaya <- function(atlas_filename, roi_numba) {
f <- function(nii_filename) {
return(paste0(bd_atlases,"/", nii_filename))
}
nii_atlas <- RNifti::readNifti( f(atlas_filename) )
# view(nii_atlas)
nii_atlas_thr <- ifelse(nii_atlas == roi_numba, 1, 0)
RNifti::writeNifti(nii_atlas_thr, template = nii_atlas, paste0(bd_atlases,"/tmp_ROI.nii.gz"))
papaya(
c(f("Dummy.nii.gz"), f("MNI152_T1_2mm_brain.nii.gz"), f("tmp_ROI.nii.gz")),
options = list(
papayaOptions(alpha = 1, lut = "Grayscale"),
papayaOptions(alpha = 0.5, lut = "Grayscale"),
papayaOptions(alpha = 0.5, lut = "Red Overlay", min = 0, max = 5)
),
interpolation = FALSE,
orthogonal = FALSE
)
}
prepare_papaya(atlas_filename, roi_numba = 6)
global_reord = FALSE
global_fontsize = 7
global_side_mm = 100
roi_numba = 6
# fmri one roi
tril_fmri <- plot_tril(
mean_RDMs_fmri[roi_numba,]$data, model = paste0("fmri_roi",roi_numba),
fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)
# emotion
tril_emotion <- plot_tril(
mean_RDMs_rats$rdm_emotion, model = "emotion",
fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)
# arousal
tril_arousal <- plot_tril(
mean_RDMs_rats$rdm_arousal, model = "arousal",
fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)
tril_fmri + tril_emotion
# We save the following into /data00/leonardo/RSA/analyses/RSA_ROI_APP
# RSA
# RSA_mean
# atlas_filename
# mean_RDMs_fmri
# mean_RDMs_rats
create_filename_results <- function(atlas_filename, copes_type, remove_neutrals) {
atlas_root <- str_replace(atlas_filename, ".nii.gz", "")
if(str_detect(copes_type,"one_ev_per_movie")) {EV_abbr <- "_OEPM"}
if(str_detect(copes_type,"emotion_predictors")) {EV_abbr <- "_EP"}
if(str_detect(copes_type,"emotion_high_low_predictors")) {EV_abbr <- "_EHLP"}
option_MN <- ifelse(str_detect(copes_type,"_minus_neutral"),"_MN","")
option_NN <- ifelse(remove_neutrals == "YES", "_RN", "")
results_filename <- paste0(atlas_root, EV_abbr, option_MN, option_NN)
return(paste0(bd_results,"/",results_filename,".RData"))
}
results_file <- create_filename_results(atlas_filename, copes_type, remove_neutrals)
results_file
## [1] "/data00/leonardo/RSA/analyses/RSA_ROI_APP/results_RSA_ROI/Yeo17_EP_MN.RData"
# save the results
save(RSA, RSA_mean, atlas_filename, mean_RDMs_fmri, mean_RDMs_rats, df_path_copes, subs, file = results_file)
# load them again
# load(results_file)
# load("/data00/leonardo/RSA/analyses/RSA_ROI_APP/results_RSA_ROI/Yeo17_OEPM_MN_RN.RData")